library(effects)
library(jtools)
library(stargazer)
library(dotwhisker)
library(ggplot2)
library(weights)
library(AER)
library(gridExtra)
library(ggpubr)
library(formattable)
library(interactions)
library(reshape)
library(dplyr)
library(tidyr)

MOC <- function (data, its, model, prog.int=100, rsq=TRUE) {
  ## Method of Composition (Caughey and Warshaw, 2018; Treier and Jackman, 2008)
  tildeB <- as.data.frame(matrix(nrow=length(its), ncol=length(coef(model))))
  R2 <- data.frame(rsq=rep(NA, length(its)), adjrsq=rep(NA, length(its)))
  for (s in seq_along(its)) {
    if (!s %% prog.int) print(s)
    it <- its[s]
    ## (1) Sample from p(x)
    data_s <- subset(data, It == it)
    #names(data_s) <- gsub("\\_s$", "\\_mean", names(data_s)) 
    ## (2) Sample from p(B|x,y):
    ##     (a) Estimate B_s and Cov(B_s) conditional on x_s.
    mod_s <- update(model, data=data_s)
    hatB_s <- coef(mod_s)
    hatV_s <- vcov(mod_s)
    ##     (b) Sample \tilde{B_s} from MV(\hat{B_s}, \hat{Cov(B_s)}).
    tildeB[s, ] <- MASS::mvrnorm(n = 1, mu = hatB_s, Sigma = hatV_s)
    if (rsq) {
      R2$rsq[s] <- summary(mod_s)$r.squared["rsq"]
      R2$adjrsq[s] <- summary(mod_s)$r.squared["adjrsq"]
    }
  }
  names(tildeB) <- names(coef(model))
  if (rsq) {
    tildeB$rsq <- R2$rsq
    tildeB$adjrsq <- R2$adjrsq
  }
  return(tildeB)
}

##### ANALYSIS

dataend2 <- read.csv("replicationdatapoqfinal.csv")

# Figure 3

theme_set(theme_bw())
all <- ggplot(dataend2, aes(x=mrp_mean, y=allidealMCMCmean)) + 
  geom_point(col="gray") +
  scale_y_continuous(expand=c(0,0.05), limits=c(-2.5,2.5)) +
  xlim(c(-1.5,1.5))+
  geom_smooth(method="lm",se=T, col="blue",fullrange=T) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") +
  ggtitle("All Votes") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  coord_cartesian(xlim=c(-1,1), ylim=c(-2,2))
econ <- ggplot(dataend2, aes(x=econmrpMCMCmean, y=econidealMCMCmean)) + 
  geom_point(col="gray") +
  scale_y_continuous(expand=c(0,0.05), limits=c(-2.5,2.5)) +
  xlim(c(-1.5,1.5))+
  geom_smooth(method="lm",se=T, col="blue",fullrange=T) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") +
  ggtitle("Economic/Spending") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(xlim=c(-1,1), ylim=c(-2,2))
env <- ggplot(dataend2, aes(x=envmrpMCMCmean, y=envidealMCMCmean)) + 
  geom_point(col="gray") +
  scale_y_continuous(expand=c(0,0.05), limits=c(-3,3)) +
  xlim(c(-1.5,1.5))+
  geom_smooth(method="lm",se=T, col="blue",fullrange=T) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") +
  ggtitle("Environment") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(xlim=c(-1,1), ylim=c(-2,2))
crim <- ggplot(dataend2, aes(x=crimmrpMCMCmean, y=crimidealMCMCmean)) + 
  geom_point(col="gray") +
  scale_y_continuous(expand=c(0,0.05), limits=c(-2.5,2.5)) +
  xlim(c(-1.5,1.5))+
  geom_smooth(method="lm",se=T, col="blue",fullrange=T) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") +
  ggtitle("Crime") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_cartesian(xlim=c(-1,1), ylim=c(-2,2))
tiff("Pomirchy-Figure 3",width=1400,res=250,height=1000)
ggarrange(all, econ,crim,env,nrow=2,ncol=2,
                    widths=c(1,1))
dev.off()

# Figure 4

ARallreg <- lm(allidealMCMCmean~mrp_mean,data=dataend2[dataend2$StateAbbr=="AR",])
theme_set(theme_bw())
ARplot <- ggplot(dataend2[dataend2$StateAbbr=="AR",], aes(x=mrp_mean, y=allidealMCMCmean,color=as.factor(1 - Republican),
                                                              shape=as.factor(Republican+1))) + 
  geom_point() +
  geom_abline(intercept=ARallreg$coefficients["(Intercept)"],slope=ARallreg$coefficients["mrp_mean"],
              col="gray",lwd=1) + 
  xlim(c(-1,1))+
  ylim(c(-2,2))+
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") + 
  ggtitle("Arkansas") +
  scale_color_discrete(name = "",
                       breaks=c("1","0"),
                       labels=c("Democrat","Republican")) + 
  scale_shape_discrete(name  ="",
                       breaks=c("1","2"),
                       labels=c("Democrat", "Republican")) +
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(legend.position="bottom")
ILallreg <- lm(allidealMCMCmean~mrp_mean,data=dataend2[dataend2$StateAbbr=="IL",])
ILplot <- ggplot(dataend2[dataend2$StateAbbr=="IL",], aes(x=mrp_mean, y=allidealMCMCmean,color=as.factor(1 - Republican),
                                                              shape=as.factor(Republican+1))) + 
  geom_point() +
  xlim(c(-1,1))+
  ylim(c(-2,2))+
  geom_abline(intercept=ILallreg$coefficients["(Intercept)"],slope=ILallreg$coefficients["mrp_mean"],
              col="gray",lwd=1) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") + 
  ggtitle("Illinois") +
  scale_color_discrete(name = "",
                       breaks=c("1","0"),
                       labels=c("Democrat","Republican")) + 
  scale_shape_discrete(name  ="",
                       breaks=c("1","2"),
                       labels=c("Democrat", "Republican")) +
  theme(plot.title = element_text(hjust = 0.5))
FLallreg <- lm(allidealMCMCmean~mrp_mean,data=dataend2[dataend2$StateAbbr=="FL",])
FLplot <- ggplot(dataend2[dataend2$StateAbbr=="FL",], aes(x=mrp_mean, y=allidealMCMCmean,color=as.factor(1 - Republican),
                                                              shape=as.factor(Republican+1))) + 
  geom_point() +
  xlim(c(-1,1))+
  ylim(c(-2,2))+
  geom_abline(intercept=FLallreg$coefficients["(Intercept)"],slope=FLallreg$coefficients["mrp_mean"],
              col="gray",lwd=1) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point") + 
  scale_color_discrete(name = "",
                       breaks=c("1","0"),
                       labels=c("Democrat","Republican")) + 
  scale_shape_discrete(name  ="",
                       breaks=c("1","2"),
                       labels=c("Democrat", "Republican")) +
  ggtitle("Florida") +
  theme(plot.title = element_text(hjust = 0.5))
TXallreg <- lm(allidealMCMCmean~mrp_mean,data=dataend2[dataend2$StateAbbr=="TX",])
TXplot <- ggplot(dataend2[dataend2$StateAbbr=="TX",], aes(x=mrp_mean, y=allidealMCMCmean,color=as.factor(1 - Republican),
                                                              shape=as.factor(Republican+1))) + 
  geom_point()+
  xlim(c(-1,1))+
  ylim(c(-2,2))+
  geom_abline(intercept=TXallreg$coefficients["(Intercept)"],slope=TXallreg$coefficients["mrp_mean"],
              col="gray",lwd=1) + 
  labs(y="Legislator Ideal Point", 
       x="Constituency Ideal Point")+ 
  scale_color_discrete(name = "",
                       breaks=c("1","0"),
                       labels=c("Democrat","Republican")) + 
  scale_shape_discrete(name  ="",
                       breaks=c("1","2"),
                       labels=c("Democrat", "Republican")) +
  ggtitle("Texas") +
  theme(plot.title = element_text(hjust = 0.5))

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
mylegend <- g_legend(ARplot)
tiff("Pomirchy-Figure 4",width=1400,res=250,height=1000)
grid.arrange(arrangeGrob(ARplot + theme(legend.position="none"),
                         ILplot + theme(legend.position="none"),
                         TXplot + theme(legend.position="none"),
                         FLplot + theme(legend.position="none"),
                         nrow=2,ncol=2),
             mylegend, nrow=2,heights=c(10, 1))
dev.off()

# Figure 5

allidealMCMC <- read.csv("allidealMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
econidealMCMC <- read.csv("econidealMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
envidealMCMC <- read.csv("envidealMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
crimidealMCMC <- read.csv("crimidealMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
allmrpMCMC <- read.csv("allmrpMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
econmrpMCMC <- read.csv("econmrpMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
envmrpMCMC <- read.csv("envmrpMCMC.csv") %>% select(num_range(prefix="X",range=1:500))
crimmrpMCMC <- read.csv("crimmrpMCMC.csv") %>% select(num_range(prefix="X",range=1:500))

alldataideal <- data.frame(cbind(dataend2,allidealMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
alldatamrp <- data.frame(cbind(dataend2,allmrpMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
    melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
alldata <- data.frame(cbind(alldataideal,alldatamrp[,"value"]))
colnames(alldata) <- c("treattwo","Republican","StateAbbr","Session",
                       "It","ideal","opinion")

econdataideal <- data.frame(cbind(dataend2,econidealMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
econdatamrp <- data.frame(cbind(dataend2,econmrpMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
econdata <- data.frame(cbind(econdataideal,econdatamrp[,"value"]))
colnames(econdata) <- c("treattwo","Republican","StateAbbr","Session",
                       "It","ideal","opinion")

envdataideal <- data.frame(cbind(dataend2,envidealMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
envdatamrp <- data.frame(cbind(dataend2,envmrpMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
envdata <- data.frame(cbind(envdataideal,envdatamrp[,"value"]))
colnames(envdata) <- c("treattwo","Republican","StateAbbr","Session",
                        "It","ideal","opinion")

crimdataideal <- data.frame(cbind(dataend2,crimidealMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
crimdatamrp <- data.frame(cbind(dataend2,crimmrpMCMC)) %>% select(num_range(prefix="X",range=1:500),treattwo,Republican,StateAbbr,Session) %>%
  melt(id.vars=c("treattwo","Republican","StateAbbr","Session"))
crimdata <- data.frame(cbind(crimdataideal,crimdatamrp[,"value"]))
colnames(crimdata) <- c("treattwo","Republican","StateAbbr","Session",
                       "It","ideal","opinion")

alldata$It <- as.numeric(gsub("X","",alldata$It)) 
econdata$It <- as.numeric(gsub("X","",econdata$It)) 
envdata$It <- as.numeric(gsub("X","",envdata$It)) 
crimdata$It <- as.numeric(gsub("X","",crimdata$It)) 

its <- sort(sample(unique(alldata$It), min(500, length(unique(alldata$It)))))
ElectoralProx <- lm(ideal~opinion*treattwo+Republican+StateAbbr+as.factor(Session),data=alldata[alldata$It==1,])
ElectoralProxEcon <- lm(ideal~opinion*treattwo+Republican+StateAbbr+as.factor(Session),data=econdata[econdata$It==1,])
ElectoralProxEnv <- lm(ideal~opinion*treattwo+Republican+StateAbbr+as.factor(Session),data=envdata[envdata$It==1,])
ElectoralProxCrim <- lm(ideal~opinion*treattwo+Republican+StateAbbr+as.factor(Session),data=crimdata[crimdata$It==1,])
set.seed(401)
moc1 <- MOC(data=alldata,its=its,model=ElectoralProx,rsq=FALSE)
moc2 <- MOC(data=econdata,its=its,model=ElectoralProxEcon,rsq=FALSE)
moc3 <- MOC(data=envdata,its=its,model=ElectoralProxEnv,rsq=FALSE)
moc4 <- MOC(data=crimdata,its=its,model=ElectoralProxCrim,rsq=FALSE)

moc.est <- list(apply(moc1,2,mean,na.rm=T),apply(moc2,2,mean,na.rm=T),
                apply(moc3,2,mean,na.rm=T),apply(moc4,2,mean,na.rm=T))
moc.sd <- list(apply(moc1,2,sd,na.rm=T),apply(moc2,2,sd,na.rm=T),
               apply(moc3,2,sd,na.rm=T),apply(moc4,2,sd,na.rm=T))

term <- c("Four Years", "Two Years", 
          "Difference in Effects")
estimate <- c(moc.est[[1]]["opinion"],
              moc.est[[1]]["opinion"]+moc.est[[1]]["opinion:treattwo"],
              moc.est[[1]]["opinion:treattwo"])
std.error <- c(moc.sd[[1]]["opinion"],
               sqrt(moc.sd[[1]]["opinion"]^2+
                      moc.sd[[1]]["treattwo"]^2+2*cov(moc1$opinion,moc1$treattwo)),
               moc.sd[[1]]["opinion:treattwo"])
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
yvalues <- rev(c(0.75,1.75,2.75))
model <- rep("All Votes",3)
dfplot <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[2]]["opinion"],
              moc.est[[2]]["opinion"]+moc.est[[2]]["opinion:treattwo"],
              moc.est[[2]]["opinion:treattwo"])
std.error <- c(moc.sd[[2]]["opinion"],
               sqrt(moc.sd[[2]]["opinion"]^2+
                      moc.sd[[2]]["treattwo"]^2+2*cov(moc1$opinion,moc1$treattwo)),
               moc.sd[[2]]["opinion:treattwo"])
model <- rep("Economic", 3)
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
yvalues <- rev(c(0.925,1.925,2.925))
dfplot2 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[3]]["opinion"],
              moc.est[[3]]["opinion"]+moc.est[[3]]["opinion:treattwo"],
              moc.est[[3]]["opinion:treattwo"])
std.error <- c(moc.sd[[3]]["opinion"],
               sqrt(moc.sd[[3]]["opinion"]^2+
                      moc.sd[[3]]["treattwo"]^2+2*cov(moc1$opinion,moc1$treattwo)),
               moc.sd[[3]]["opinion:treattwo"])
model <- rep("Environment", 3)
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
yvalues <- rev(c(1.1,2.1,3.1))
dfplot3 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[4]]["opinion"],
              moc.est[[4]]["opinion"]+moc.est[[4]]["opinion:treattwo"],
              moc.est[[4]]["opinion:treattwo"])
std.error <- c(moc.sd[[4]]["opinion"],
               sqrt(moc.sd[[4]]["opinion"]^2+
                      moc.sd[[4]]["treattwo"]^2+2*cov(moc1$opinion,moc1$treattwo)),
               moc.sd[[4]]["opinion:treattwo"])
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
model <- rep("Crime", 3)
yvalues <- rev(c(1.275,2.275,3.275))
dfplot4 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
dfplot5 <- rbind(dfplot,dfplot2,dfplot3,dfplot4)
dfplot5$model <- as.factor(dfplot5$model)
tiff("Pomirchy-Figure 5",width=1400,res=250,height=1000)
theme_set(theme_bw())
dwplot(dfplot5,dodge_size=0.7,
       model_order = c("All Votes","Economic","Environment","Crime"),
       dot_args = list(aes(shape = model))) +
  geom_vline(aes(x=0), lty=2, xintercept=0) +
  ylab('Term Length') + 
  xlab('Coefficient') +
  xlim(c(-2,2)) + 
  scale_y_discrete(name = 'Variable', 
                   labels = c('Difference \n in Effects',
                              'Two Years',
                              'Four Years')) +
  geom_text(aes(x=estimate,y=yvalues,label=pvalues),size=2,nudge_y=0.075)+
  theme(plot.title = element_text(face="bold",hjust=0.5),
        legend.direction = "horizontal",
        legend.position = "bottom")
dev.off()

# Figure 6

ElectoralProxAR <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=alldata[alldata$It==1 & alldata$StateAbbr=="AR",])
ElectoralProxFL <- lm(ideal~opinion*treattwo+Republican,data=alldata[alldata$It==1 & alldata$StateAbbr=="FL",])
ElectoralProxIL <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=alldata[alldata$It==1 & alldata$StateAbbr=="IL",])
ElectoralProxTX <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=alldata[alldata$It==1 & alldata$StateAbbr=="TX",])
ElectoralProxEconAR <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=econdata[econdata$It==1 & econdata$StateAbbr=="AR",])
ElectoralProxEconFL <- lm(ideal~opinion*treattwo+Republican,data=econdata[econdata$It==1 & econdata$StateAbbr=="FL",])
ElectoralProxEconIL <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=econdata[econdata$It==1 & econdata$StateAbbr=="IL",])
ElectoralProxEconTX <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=econdata[econdata$It==1 & econdata$StateAbbr=="TX",])
ElectoralProxEnvAR <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=envdata[envdata$It==1 & envdata$StateAbbr=="AR",])
ElectoralProxEnvFL <- lm(ideal~opinion*treattwo+Republican,data=envdata[envdata$It==1 & envdata$StateAbbr=="FL",])
ElectoralProxEnvIL <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=envdata[envdata$It==1 & envdata$StateAbbr=="IL",])
ElectoralProxEnvTX <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=envdata[envdata$It==1 & envdata$StateAbbr=="TX",])
ElectoralProxCrimAR <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=crimdata[crimdata$It==1 & crimdata$StateAbbr=="AR",])
ElectoralProxCrimFL <- lm(ideal~opinion*treattwo+Republican,data=crimdata[crimdata$It==1 & crimdata$StateAbbr=="FL",])
ElectoralProxCrimIL <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=crimdata[crimdata$It==1 & crimdata$StateAbbr=="IL",])
ElectoralProxCrimTX <- lm(ideal~opinion*treattwo+Republican+as.factor(Session),data=crimdata[crimdata$It==1 & crimdata$StateAbbr=="TX",])

its <- sort(sample(unique(alldata$It), min(500, length(unique(alldata$It)))))

set.seed(401)
moc1 <- MOC(data=subset(alldata,alldata$StateAbbr=="AR"),its=its,model=ElectoralProxAR,rsq=FALSE)
moc2 <- MOC(data=subset(alldata,alldata$StateAbbr=="FL"),its=its,model=ElectoralProxFL,rsq=FALSE)
moc3 <- MOC(data=subset(alldata,alldata$StateAbbr=="IL"),its=its,model=ElectoralProxIL,rsq=FALSE)
moc4 <- MOC(data=subset(alldata,alldata$StateAbbr=="TX"),its=its,model=ElectoralProxTX,rsq=FALSE)
moc5 <- MOC(data=subset(econdata,econdata$StateAbbr=="AR"),its=its,model=ElectoralProxEconAR,rsq=FALSE)
moc6 <- MOC(data=subset(econdata,econdata$StateAbbr=="FL"),its=its,model=ElectoralProxEconFL,rsq=FALSE)
moc7 <- MOC(data=subset(econdata,econdata$StateAbbr=="IL"),its=its,model=ElectoralProxEconIL,rsq=FALSE)
moc8 <- MOC(data=subset(econdata,econdata$StateAbbr=="TX"),its=its,model=ElectoralProxEconTX,rsq=FALSE)
moc9 <- MOC(data=subset(envdata,envdata$StateAbbr=="AR"),its=its,model=ElectoralProxEnvAR,rsq=FALSE)
moc10 <- MOC(data=subset(envdata,envdata$StateAbbr=="FL"),its=its,model=ElectoralProxEnvFL,rsq=FALSE)
moc11 <- MOC(data=subset(envdata,envdata$StateAbbr=="IL"),its=its,model=ElectoralProxEnvIL,rsq=FALSE)
moc12 <- MOC(data=subset(envdata,envdata$StateAbbr=="TX"),its=its,model=ElectoralProxEnvTX,rsq=FALSE)
moc13 <- MOC(data=subset(crimdata,crimdata$StateAbbr=="AR"),its=its,model=ElectoralProxCrimAR,rsq=FALSE)
moc14 <- MOC(data=subset(crimdata,crimdata$StateAbbr=="FL"),its=its,model=ElectoralProxCrimFL,rsq=FALSE)
moc15 <- MOC(data=subset(crimdata,crimdata$StateAbbr=="IL"),its=its,model=ElectoralProxCrimIL,rsq=FALSE)
moc16 <- MOC(data=subset(crimdata,crimdata$StateAbbr=="TX"),its=its,model=ElectoralProxCrimTX,rsq=FALSE)

moc.est <- list(apply(moc1,2,mean,na.rm=T),apply(moc2,2,mean,na.rm=T),
                apply(moc3,2,mean,na.rm=T),apply(moc4,2,mean,na.rm=T),
                apply(moc5,2,mean,na.rm=T),apply(moc6,2,mean,na.rm=T),
                apply(moc7,2,mean,na.rm=T),apply(moc8,2,mean,na.rm=T),
                apply(moc9,2,mean,na.rm=T),apply(moc10,2,mean,na.rm=T),
                apply(moc11,2,mean,na.rm=T),apply(moc12,2,mean,na.rm=T),
                apply(moc13,2,mean,na.rm=T),apply(moc14,2,mean,na.rm=T),
                apply(moc15,2,mean,na.rm=T),apply(moc16,2,mean,na.rm=T))
moc.sd <- list(apply(moc1,2,sd,na.rm=T),apply(moc2,2,sd,na.rm=T),
               apply(moc3,2,sd,na.rm=T),apply(moc4,2,sd,na.rm=T),
               apply(moc5,2,sd,na.rm=T),apply(moc6,2,sd,na.rm=T),
               apply(moc7,2,sd,na.rm=T),apply(moc8,2,sd,na.rm=T),
               apply(moc9,2,sd,na.rm=T),apply(moc10,2,sd,na.rm=T),
               apply(moc11,2,sd,na.rm=T),apply(moc12,2,sd,na.rm=T),
               apply(moc13,2,sd,na.rm=T),apply(moc14,2,sd,na.rm=T),
               apply(moc15,2,sd,na.rm=T),apply(moc16,2,sd,na.rm=T))

term <- c("Arkansas", "Florida", 
          "Illinois", "Texas")
estimate <- c(moc.est[[1]]["opinion:treattwo"],
              moc.est[[2]]["opinion:treattwo"],
              moc.est[[3]]["opinion:treattwo"],
              moc.est[[4]]["opinion:treattwo"])
std.error <- c(moc.sd[[1]]["opinion:treattwo"],
               moc.sd[[2]]["opinion:treattwo"],
               moc.sd[[3]]["opinion:treattwo"],
               moc.sd[[4]]["opinion:treattwo"])
model <- rep("All Votes",4)
yvalues <- rev(c(0.75,1.75,2.75,3.75))
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
dfplot <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[5]]["opinion:treattwo"],
              moc.est[[6]]["opinion:treattwo"],
              moc.est[[7]]["opinion:treattwo"],
              moc.est[[8]]["opinion:treattwo"])
std.error <- c(moc.sd[[5]]["opinion:treattwo"],
               moc.sd[[6]]["opinion:treattwo"],
               moc.sd[[7]]["opinion:treattwo"],
               moc.sd[[8]]["opinion:treattwo"])
model <- rep("Economic", 4)
yvalues <- rev(c(0.925,1.925,2.925,3.925))
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
dfplot2 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[9]]["opinion:treattwo"],
              moc.est[[10]]["opinion:treattwo"],
              moc.est[[11]]["opinion:treattwo"],
              moc.est[[12]]["opinion:treattwo"])
std.error <- c(moc.sd[[9]]["opinion:treattwo"],
               moc.sd[[10]]["opinion:treattwo"],
               moc.sd[[11]]["opinion:treattwo"],
               moc.sd[[12]]["opinion:treattwo"])
model <- rep("Environment", 4)
yvalues <- rev(c(1.1,2.1,3.1,4.1))
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
dfplot3 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
estimate <- c(moc.est[[13]]["opinion:treattwo"],
              moc.est[[14]]["opinion:treattwo"],
              moc.est[[15]]["opinion:treattwo"],
              moc.est[[16]]["opinion:treattwo"])
std.error <- c(moc.sd[[13]]["opinion:treattwo"],
               moc.sd[[14]]["opinion:treattwo"],
               moc.sd[[15]]["opinion:treattwo"],
               moc.sd[[16]]["opinion:treattwo"])
model <- rep("Crime", 4)
yvalues <- rev(c(1.275,2.275,3.275,4.275))
pvalues <- round(2*pnorm(-abs(estimate/std.error)),3)
pvalues <- ifelse(pvalues==0,paste("p < 0.001"),paste("p=",pvalues,sep=""))
dfplot4 <- data.frame(term, estimate, std.error, model,pvalues,yvalues)
dfplot5 <- rbind(dfplot,dfplot2,dfplot3,dfplot4)
dfplot5$model <- as.factor(dfplot5$model)
tiff("Pomirchy-Figure 6",width=1400,res=250,height=1000)
theme_set(theme_bw())
dwplot(dfplot5,dodge_size=0.7,
       dot_args = list(aes(shape = model))) +
  geom_vline(aes(x=0), lty=2, xintercept=0) +
  ylab('Results') + 
  xlab('Coefficient') + 
  xlim(c(-4,4))+
  scale_y_discrete(name = 'State', 
                   labels = c('Texas',
                              'Illinois',
                              'Florida',
                              'Arkansas')) +
  geom_text(aes(x=estimate,y=yvalues,label=pvalues),size=2,nudge_y=0.085)+
  theme(plot.title = element_text(face="bold",hjust=0.5),
        legend.direction = "horizontal",
        legend.position = "bottom")
dev.off()

# Figure 7

dataIL <- dataend2[dataend2$StateAbbr=="IL" & dataend2$Session==2012,]
dataIL$condition <- as.factor(ifelse(dataIL$treattwo==1,"Two Years",
                                     "Four Years"))

gaysuppopinion <- read.csv("IllinoisestimatesGayMarriageMRP.csv")
colnames(gaysuppopinion)[which(colnames(gaysuppopinion)=="X")] <- "DistrictNumber"

gaysuppdata <- merge(dataIL,gaysuppopinion,by="DistrictNumber")
gaysuppdata$vote <- ifelse(gaysuppdata$GaySupp==1,1,
                           ifelse(gaysuppdata$GaySupp==6,0,NA))
gaysuppdata$overall <- gaysuppdata$mean
gaysuppdata$majsupport <- ifelse(gaysuppdata$overall>=0.5,1,0)
gaysuppdata$resp <- gaysuppdata$vote==gaysuppdata$majsupport

GaySupp <- lm(vote~overall*condition,data=gaysuppdata)

penreformopinion <- read.csv("IllinoisestimatesPenReformMRP.csv")
colnames(penreformopinion)[which(colnames(penreformopinion)=="X")] <- "DistrictNumber"

penreformdata <- merge(dataIL,penreformopinion,by="DistrictNumber")
penreformdata$vote <- ifelse(penreformdata$PensionReform==1,1,
                             ifelse(penreformdata$PensionReform==6,0,NA))
penreformdata$overall <- penreformdata$mean
penreformdata$majsupport <- ifelse(penreformdata$overall>=0.5,1,0)
penreformdata$resp <- penreformdata$vote==penreformdata$majsupport

PenReform <- lm(vote~overall*condition,data=penreformdata)

medmjopinion <- read.csv("IllinoisestimatesMedMJMRP.csv")
colnames(medmjopinion)[which(colnames(medmjopinion)=="X")] <- "DistrictNumber"

medmjdata <- merge(dataIL,medmjopinion,by="DistrictNumber")
medmjdata$vote <- ifelse(medmjdata$MedMJ==1,1,
                         ifelse(medmjdata$MedMJ==6,0,NA))
medmjdata$overall <- medmjdata$mean
medmjdata$majsupport <- ifelse(medmjdata$overall>=0.5,1,0)
medmjdata$resp <- medmjdata$vote==medmjdata$majsupport

MedMJ <- lm(vote~overall*condition,data=medmjdata)

conccarryopinion <- read.csv("IllinoisestimatesConcealCarryMRP.csv")
colnames(conccarryopinion)[which(colnames(conccarryopinion)=="X")] <- "DistrictNumber"

conccarrydata <- merge(dataIL,conccarryopinion,by="DistrictNumber")
conccarrydata$vote <- ifelse(conccarrydata$ConcealCarry==1,1,
                             ifelse(conccarrydata$ConcealCarry==6,0,NA))
conccarrydata$overall <- conccarrydata$mean
conccarrydata$majsupport <- ifelse(conccarrydata$overall>=0.5,1,0)
conccarrydata$resp <- conccarrydata$vote==conccarrydata$majsupport

ConcCarry <- lm(vote~overall*condition,data=conccarrydata)

gamblingopinion <- read.csv("IllinoisestimatesGambleMRP.csv")
colnames(gamblingopinion)[which(colnames(gamblingopinion)=="X")] <- "DistrictNumber"

gamblingdata <- merge(dataIL,gamblingopinion,by="DistrictNumber")
gamblingdata$vote <- ifelse(gamblingdata$Gamble==1,1,
                            ifelse(gamblingdata$Gamble==6,0,NA))
gamblingdata$overall <- gamblingdata$mean
gamblingdata$majsupport <- ifelse(gamblingdata$overall>=0.5,1,0)
gamblingdata$resp <- gamblingdata$vote==gamblingdata$majsupport

Gambling <- lm(vote~overall*condition,data=gamblingdata)

checksopinion <- read.csv("IllinoisestimatesBackgroundChecksMRP.csv")
colnames(checksopinion)[which(colnames(checksopinion)=="X")] <- "DistrictNumber"

checksdata <- merge(dataIL,checksopinion,by="DistrictNumber")
checksdata$vote <- ifelse(checksdata$Background==1,1,
                          ifelse(checksdata$Background==6,0,NA))
checksdata$overall <- checksdata$mean
checksdata$majsupport <- ifelse(checksdata$overall>=0.5,1,0)
checksdata$resp <- checksdata$vote==checksdata$majsupport

Checks <- lm(vote~overall*condition,data=checksdata)

cola25Kopinion <- read.csv("IllinoisestimatesCOLA25KMRP.csv")
colnames(cola25Kopinion)[which(colnames(cola25Kopinion)=="X")] <- "DistrictNumber"

cola25Kdata <- merge(dataIL,cola25Kopinion,by="DistrictNumber")
cola25Kdata$vote <- ifelse(cola25Kdata$COLA25K==1,1,
                           ifelse(cola25Kdata$COLA25K==6,0,NA))
cola25Kdata$overall <- cola25Kdata$mean
cola25Kdata$majsupport <- ifelse(cola25Kdata$overall>=0.5,1,0)
cola25Kdata$resp <- cola25Kdata$vote==cola25Kdata$majsupport

cola25K <- lm(vote~overall*condition,data=cola25Kdata)

ammoopinion <- read.csv("IllinoisestimatesAmmoMRP.csv")
colnames(ammoopinion)[which(colnames(ammoopinion)=="X")] <- "DistrictNumber"

ammodata <- merge(dataIL,ammoopinion,by="DistrictNumber")
ammodata$vote <- ifelse(ammodata$Ammo==1,1,
                        ifelse(ammodata$Ammo==6,0,NA))
ammodata$overall <- ammodata$mean
ammodata$majsupport <- ifelse(ammodata$overall>=0.5,1,0)
ammodata$resp <- ammodata$vote==ammodata$majsupport

ammo <- lm(vote~overall*condition,data=ammodata)

fulldata <- rbind(ammodata,checksdata,conccarrydata,gamblingdata,
                  gaysuppdata,medmjdata,penreformdata,cola25Kdata)
fulldata$issue <- c(rep(0,nrow(ammodata)),rep(1,nrow(checksdata)),
                    rep(2,nrow(conccarrydata)),rep(3,nrow(gamblingdata)),
                    rep(4,nrow(gaysuppdata)),rep(5,nrow(medmjdata)),
                    rep(6,nrow(penreformdata)),rep(7,nrow(cola25Kdata)))
fullreg <- lm(vote~overall*condition,data=fulldata)

ammoplot <- interact_plot(ammo, pred = "overall", modx = "condition", 
                          x.label="Constituency Opinion", y.label="Vote in Favor",
                          plot.points=T,legend.main="Condition") + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Caps on Ammunition") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom")
checksplot <- interact_plot(Checks, pred = "overall", modx = "condition", 
                            x.label="Constituency Opinion", y.label="Vote in Favor",
                            plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Background Checks") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
gambleplot <- interact_plot(Gambling, pred = "overall", modx = "condition", 
                            x.label="Constituency Opinion", y.label="Vote in Favor",
                            plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Expand Gambling") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
medmjplot <- interact_plot(MedMJ, pred = "overall", modx = "condition", 
                           x.label="Constituency Opinion", y.label="Vote in Favor",
                           plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Med. Marijuana") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
penreformplot <- interact_plot(PenReform, pred = "overall", modx = "condition", 
                               x.label="Constituency Opinion", y.label="Vote in Favor",
                               plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Pension Reform") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
cola25Kplot <- interact_plot(cola25K, pred = "overall", modx = "condition", 
                             x.label="Constituency Opinion", y.label="Vote in Favor",
                             plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("COLA Increases") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
gaysuppplot <- interact_plot(GaySupp, pred = "overall", modx = "condition", 
                             x.label="Constituency Opinion", y.label="Vote in Favor",
                             plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Same-Sex Marriage") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
conccarryplot <- interact_plot(ConcCarry, pred = "overall", modx = "condition", 
                               x.label="Constituency Opinion", y.label="Vote in Favor",
                               plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("Expand Conc. Carry") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
pooledplot <- interact_plot(fullreg, pred = "overall", modx = "condition", 
                            x.label="Constituency Opinion", y.label="Vote in Favor",
                            plot.points=T) + 
  ylim(c(0,1)) + 
  xlim(c(0,1)) + 
  ggtitle("All Votes") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position = "none")
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
mylegend <- g_legend(ammoplot)
ammoplot <- ammoplot + theme(legend.position="none")
tiff("Pomirchy-Figure 7",width=2000,res=250,height=2000)
grid.arrange(arrangeGrob(pooledplot,ammoplot,checksplot,
                         conccarryplot,gambleplot,gaysuppplot,
                         medmjplot,penreformplot,cola25Kplot,
                         nrow=3,ncol=3),
             mylegend,nrow=2,heights=c(11,1))
dev.off()

# Table 3

gaysuppopinion <- read.csv("IllinoisestimatesGayMarriageMRP.csv")
colnames(gaysuppopinion)[which(colnames(gaysuppopinion)=="X")] <- "DistrictNumber"

dataIL <- dataend2[dataend2$StateAbbr=="IL" & dataend2$Session==2012,]
gaysuppdata <- merge(dataIL,gaysuppopinion,by="DistrictNumber")
gaysuppdata$vote <- ifelse(gaysuppdata$GaySupp==1,1,
                           ifelse(gaysuppdata$GaySupp==6,0,NA))
gaysuppdata2 <- gaysuppdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(gaysuppdata2) <- c("treattwo","vote","It","opinion")
gaysuppdata2$It <- as.numeric(gsub("X","",gaysuppdata2$It))

its <- sort(sample(unique(gaysuppdata2$It), min(500, length(unique(gaysuppdata2$It)))))

gaysuppreg <- lm(vote~opinion*treattwo,data=gaysuppdata2[gaysuppdata2$It==1,])
set.seed(401)
moc1 <- MOC(data=gaysuppdata2,its=its,model=gaysuppreg,rsq=FALSE)

penreformopinion <- read.csv("IllinoisestimatesPenReformMRP.csv")
colnames(penreformopinion)[which(colnames(penreformopinion)=="X")] <- "DistrictNumber"

penreformdata <- merge(dataIL,penreformopinion,by="DistrictNumber")
penreformdata$vote <- ifelse(penreformdata$PensionReform==1,1,
                             ifelse(penreformdata$PensionReform==6,0,NA))

penreformdata2 <- penreformdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(penreformdata2) <- c("treattwo","vote","It","opinion")
penreformdata2$It <- as.numeric(gsub("X","",penreformdata2$It))

penreformreg <- lm(vote~opinion*treattwo,data=penreformdata2[penreformdata2$It==1,])
set.seed(401)
moc2 <- MOC(data=penreformdata2,its=its,model=penreformreg,rsq=FALSE)

medmjopinion <- read.csv("IllinoisestimatesMedMJMRP.csv")
colnames(medmjopinion)[which(colnames(medmjopinion)=="X")] <- "DistrictNumber"

medmjdata <- merge(dataIL,medmjopinion,by="DistrictNumber")
medmjdata$vote <- ifelse(medmjdata$MedMJ==1,1,
                         ifelse(medmjdata$MedMJ==6,0,NA))

medmjdata2 <- medmjdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(medmjdata2) <- c("treattwo","vote","It","opinion")
medmjdata2$It <- as.numeric(gsub("X","",medmjdata2$It))

medmjreg <- lm(vote~opinion*treattwo,data=medmjdata2[medmjdata2$It==1,])
set.seed(401)
moc3 <- MOC(data=medmjdata2,its=its,model=medmjreg,rsq=FALSE)

conccarryopinion <- read.csv("IllinoisestimatesConcealCarryMRP.csv")
colnames(conccarryopinion)[which(colnames(conccarryopinion)=="X")] <- "DistrictNumber"

conccarrydata <- merge(dataIL,conccarryopinion,by="DistrictNumber")
conccarrydata$vote <- ifelse(conccarrydata$ConcealCarry==1,1,
                             ifelse(conccarrydata$ConcealCarry==6,0,NA))

conccarrydata2 <- conccarrydata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(conccarrydata2) <- c("treattwo","vote","It","opinion")
conccarrydata2$It <- as.numeric(gsub("X","",conccarrydata2$It))

conccarryreg <- lm(vote~opinion*treattwo,data=conccarrydata2[conccarrydata2$It==1,])
set.seed(401)
moc4 <- MOC(data=conccarrydata2,its=its,model=conccarryreg,rsq=FALSE)

gamblingopinion <- read.csv("IllinoisestimatesGambleMRP.csv")
colnames(gamblingopinion)[which(colnames(gamblingopinion)=="X")] <- "DistrictNumber"

gamblingdata <- merge(dataIL,gamblingopinion,by="DistrictNumber")
gamblingdata$vote <- ifelse(gamblingdata$Gamble==1,1,
                            ifelse(gamblingdata$Gamble==6,0,NA))

gamblingdata2 <- gamblingdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(gamblingdata2) <- c("treattwo","vote","It","opinion")
gamblingdata2$It <- as.numeric(gsub("X","",gamblingdata2$It))

gamblingreg <- lm(vote~opinion*treattwo,data=gamblingdata2[gamblingdata2$It==1,])
set.seed(401)
moc5 <- MOC(data=gamblingdata2,its=its,model=gamblingreg,rsq=FALSE)

checksopinion <- read.csv("IllinoisestimatesBackgroundChecksMRP.csv")
colnames(checksopinion)[which(colnames(checksopinion)=="X")] <- "DistrictNumber"

checksdata <- merge(dataIL,checksopinion,by="DistrictNumber")
checksdata$vote <- ifelse(checksdata$Background==1,1,
                          ifelse(checksdata$Background==6,0,NA))

checksdata2 <- checksdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(checksdata2) <- c("treattwo","vote","It","opinion")
checksdata2$It <- as.numeric(gsub("X","",checksdata2$It))

checksreg <- lm(vote~opinion*treattwo,data=checksdata2[checksdata2$It==1,])
set.seed(401)
moc6 <- MOC(data=checksdata2,its=its,model=checksreg,rsq=FALSE)

cola25Kopinion <- read.csv("IllinoisestimatesCOLA25KMRP.csv")
colnames(cola25Kopinion)[which(colnames(cola25Kopinion)=="X")] <- "DistrictNumber"

cola25Kdata <- merge(dataIL,cola25Kopinion,by="DistrictNumber")
cola25Kdata$vote <- ifelse(cola25Kdata$COLA25K==1,1,
                           ifelse(cola25Kdata$COLA25K==6,0,NA))

cola25Kdata2 <- cola25Kdata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(cola25Kdata2) <- c("treattwo","vote","It","opinion")
cola25Kdata2$It <- as.numeric(gsub("X","",cola25Kdata2$It))

cola25Kreg <- lm(vote~opinion*treattwo,data=cola25Kdata2[cola25Kdata2$It==1,])
set.seed(401)
moc7 <- MOC(data=cola25Kdata2,its=its,model=cola25Kreg,rsq=FALSE)

ammoopinion <- read.csv("IllinoisestimatesAmmoMRP.csv")
colnames(ammoopinion)[which(colnames(ammoopinion)=="X")] <- "DistrictNumber"

ammodata <- merge(dataIL,ammoopinion,by="DistrictNumber")
ammodata$vote <- ifelse(ammodata$Ammo==1,1,
                        ifelse(ammodata$Ammo==6,0,NA))

ammodata2 <- ammodata %>% select(num_range(prefix="X",range=1:500),treattwo,vote) %>%
  melt(id.vars=c("treattwo","vote"))
colnames(ammodata2) <- c("treattwo","vote","It","opinion")
ammodata2$It <- as.numeric(gsub("X","",ammodata2$It))

ammoreg <- lm(vote~opinion*treattwo,data=ammodata2[ammodata2$It==1,])
set.seed(401)
moc8 <- MOC(data=ammodata2,its=its,model=ammoreg,rsq=FALSE)

fulldata <- rbind(ammodata2,checksdata2,conccarrydata2,gamblingdata2,
                  gaysuppdata2,medmjdata2,penreformdata2,cola25Kdata2)
fulldata$issue <- c(rep(0,nrow(ammodata2)),rep(1,nrow(checksdata2)),
                    rep(2,nrow(conccarrydata2)),rep(3,nrow(gamblingdata2)),
                    rep(4,nrow(gaysuppdata2)),rep(5,nrow(medmjdata2)),
                    rep(6,nrow(penreformdata2)),rep(7,nrow(cola25Kdata2)))
fullreg <- lm(vote~opinion*treattwo+as.factor(issue),data=fulldata[fulldata$It==1,])
set.seed(401)
moc9 <- MOC(data=fulldata,its=its,model=fullreg,rsq=FALSE)

moc.est <- list(apply(moc9,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc1,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc2,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc4,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")])
moc.est2 <- list(apply(moc3,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                 apply(moc6,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                 apply(moc7,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                 apply(moc8,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                 apply(moc5,2,mean,na.rm=T)[c("opinion","treattwo","opinion:treattwo")])
moc.sd <- list(apply(moc9,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
               apply(moc1,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
               apply(moc2,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
               apply(moc4,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")])
moc.sd2 <- list(apply(moc3,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc6,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc7,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc8,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")],
                apply(moc5,2,sd,na.rm=T)[c("opinion","treattwo","opinion:treattwo")])

stargazer(
  fullreg,gaysuppreg,penreformreg,conccarryreg,
  font.size="footnotesize",model.numbers=F,
  column.sep.width="1pt",single.row=T,type="html",out="table3.htm",
  keep=c("opinion","treattwo","opinion:treattwo"),
  covariate.labels=c("Constituency Opinion","Two-Year","Const. Opinion * Two-Year"),
  coef=moc.est,se=moc.sd,p.auto=T,report=c("vcsp"),
  title="Responsiveness on Specific Roll-Call Votes",
  dep.var.labels="Vote in Favor",
  column.labels=c("Pooled","Gay Marriage","Pension Reform","Concealed Carry"),
  omit.stat=c("F", "rsq","adj.rsq","ser")
)

## Table 4

stargazer(
  medmjreg,checksreg,cola25Kreg,ammoreg,gamblingreg,
  font.size="footnotesize",model.numbers=F,
  column.sep.width="1pt",single.row=T,type="html",out="table4.htm",
  keep=c("opinion","treattwo","opinion:treattwo"),
  covariate.labels=c("Constituency Opinion","Two-Year","Const. Opinion * Two-Year"),
  coef=moc.est2,se=moc.sd2,p.auto=T,report=c("vcsp"),
  title="Responsiveness on Specific Roll-Call Votes (cont'd)",
  dep.var.labels="Vote in Favor",
  column.labels=c("Medical Marijuana","Background Checks","COLA Increases",
                  "Ammunition Clips","Gambling"),
  omit.stat=c("F", "rsq","adj.rsq","ser")
)



